C $TITLE: 'EKSC'
C $NOFLOATCALLS
C
      SUBROUTINE EKSC (S,Z,RH,XK,IJ,EZS,ERS,EZC,ERC,EZK,ERK)
C     COMPUTE E FIELD OF SINE, COSINE, AND CONSTANT CURRENT FILAMENTS BY
C     THIN WIRE APPROXIMATION.
      REAL*8 CONX(2),XK,RKB2,RHK,SHK,SS,CS
      COMPLEX*16 CON,GZ1,GZ2,GP1,GP2,GZP1,GZP2,EZS,ERS,EZC,ERC,EZK,ERK
      COMMON /TMI/ ZPK,RKB2,IJX
      EQUIVALENCE (CONX,CON)
      DATA CONX/0.,4.771341189D0/
C**
      IJX=IJ
      ZPK=XK*Z
      RHK=XK*RH
      RKB2=RHK*RHK
      SH=.5*S
      SHK=XK*SH
      SS=DSIN(SHK)
      CS=DCOS(SHK)
      Z2=SH-Z
      Z1=-(SH+Z)
C**
      CALL GX (Z1,RH,XK,GZ1,GP1)
      CALL GX (Z2,RH,XK,GZ2,GP2)
C**
      GZP1=GP1*Z1
      GZP2=GP2*Z2
      EZS=CON*((GZ2-GZ1)*CS*XK-(GZP2+GZP1)*SS)
      EZC=-CON*((GZ2+GZ1)*SS*XK+(GZP2-GZP1)*CS)
      ERK=CON*(GP2-GP1)*RH
      CALL INTX (-SHK,SHK,RHK,IJ,CINT,SINT)
      EZK=-CON*(GZP2-GZP1+XK*XK*CMPLX(CINT,-SINT))
      GZP1=GZP1*Z1
      GZP2=GZP2*Z2
      IF (RH.LT.1.E-10) GO TO 1
      ERS=-CON*((GZP2+GZP1+GZ2+GZ1)*SS-(Z2*GZ2-Z1*GZ1)*CS*XK)/RH
      ERC=-CON*((GZP2-GZP1+GZ2-GZ1)*CS+(Z2*GZ2+Z1*GZ1)*SS*XK)/RH
C**
      RETURN
1     ERS=(0.,0.)
      ERC=(0.,0.)
      RETURN
      END
C
C
C
      SUBROUTINE EKSCX(BX,S,Z,RHX,XK,IJ,INX1,INX2,EZS,ERS,EZC,ERC,EZK,
     1ERK)
C     COMPUTE E FIELD OF SINE, COSINE, AND CONSTANT CURRENT FILAMENTS BY
C     EXTENDED THIN WIRE APPROXIMATION.
      INTEGER*4 INX1,INX2
      REAL*8 CONX(2),XK,RKB2,RHK,SHK,SS,CS
      COMPLEX*16 CON,EZS,ERS,EZC,ERC,EZK,ERK,GZ1,GZ2,GRK1,GRK2
      COMPLEX*16 GZP1,GZP2,GR1,GR2,GRP1,GRP2,GZZ1,GZZ2
      COMMON /TMI/ ZPK,RKB2,IJX
      EQUIVALENCE (CONX,CON)
      DATA CONX/0.,4.771341189D0/
C**
      IF (RHX.LT.BX) GO TO 1
      RH=RHX
      B=BX
      IRA=0
      GO TO 2
1     RH=BX
      B=RHX
      IRA=1
2     SH=.5*S
      IJX=IJ
      ZPK=XK*Z
      RHK=XK*RH
      RKB2=RHK*RHK
      SHK=XK*SH
      SS=DSIN(1.D0*SHK)
      CS=DCOS(1.D0*SHK)
      Z2=SH-Z
      Z1=-(SH+Z)
      A2=B*B
      IF (INX1.EQ.2) GO TO 3
      CALL GXX (Z1,RH,B,A2,XK,IRA,GZ1,GZP1,GR1,GRP1,GRK1,GZZ1)
      GO TO 4
3     CALL GX (Z1,RHX,XK,GZ1,GRK1)
      GZP1=GRK1*Z1
      GR1=GZ1/RHX
      GRP1=GZP1/RHX
      GRK1=GRK1*RHX
      GZZ1=(0.,0.)
4     IF (INX2.EQ.2) GO TO 5
      CALL GXX (Z2,RH,B,A2,XK,IRA,GZ2,GZP2,GR2,GRP2,GRK2,GZZ2)
      GO TO 6
5     CALL GX (Z2,RHX,XK,GZ2,GRK2)
      GZP2=GRK2*Z2
      GR2=GZ2/RHX
      GRP2=GZP2/RHX
      GRK2=GRK2*RHX
      GZZ2=(0.,0.)
6     EZS=CON*((GZ2-GZ1)*CS*XK-(GZP2+GZP1)*SS)
      EZC=-CON*((GZ2+GZ1)*SS*XK+(GZP2-GZP1)*CS)
      ERS=-CON*((Z2*GRP2+Z1*GRP1+GR2+GR1)*SS-(Z2*GR2-Z1*GR1)*CS*XK)
      ERC=-CON*((Z2*GRP2-Z1*GRP1+GR2-GR1)*CS+(Z2*GR2+Z1*GR1)*SS*XK)
      ERK=CON*(GRK2-GRK1)
      CALL INTX (-SHK,SHK,RHK,IJ,CINT,SINT)
      BK=B*XK
      BK2=BK*BK*.25
      EZK=-CON*(GZP2-GZP1+XK*XK*(1.-BK2)*CMPLX(CINT,-SINT)-BK2*(GZZ2
     1 -GZZ1))
      RETURN
      END
C
C
C
      SUBROUTINE ROM2 (A,B,SUM,DMIN)
C
C     FOR THE SOMMERFELD GROUND OPTION, ROM2 INTEGRATES OVER THE SOURCE
C     SEGMENT TO OBTAIN THE TOTAL FIELD DUE TO GROUND.  THE METHOD OF
C     VARIABLE INTERVAL WIDTH ROMBERG INTEGRATION IS USED.  THERE ARE 9
C     FIELD COMPONENTS - THE X, Y, AND Z COMPONENTS DUE TO CONSTANT,
C     SINE, AND COSINE CURRENT DISTRIBUTIONS.
C
      REAL*8 A,B,TMAG1,TMAG2,TR,TI,DMIN,Z,DZ,DZOT,S,EP,ZE,ZEND
      COMPLEX*16 SUM,G1,G2,G3,G4,G5,T00,T01,T10,T02,T11,T20
      INTEGER*4 NM,NS
      DIMENSION SUM(9),G1(9),G2(9),G3(9),G4(9),G5(9),T01(9),T10(9),
     1T20(9)
      DATA NM,NTS,NX,N/65536,4,1,9/,RX/1.E-4/
C**
      Z=A
      ZE=B
      S=B-A
      IF (S.GE.0.) GO TO 1
      WRITE(*,18)
      STOP
1     EP=S/(1.D4*NM)
      ZEND=ZE-EP
      DO 2 I=1,N
2     SUM(I)=(0.,0.)
      NS=NX
      NT=0
      CALL SFLDS (Z,G1)
3     DZ=S/NS
      IF (Z+DZ.LE.ZE) GO TO 4
      DZ=ZE-Z
      IF (DZ.LE.EP) GO TO 17
4     DZOT=DZ*.5
      CALL SFLDS (Z+DZOT,G3)
      CALL SFLDS (Z+DZ,G5)
5     TMAG1=0.
      TMAG2=0.
C
C     EVALUATE 3 POINT ROMBERG RESULT AND TEST CONVERGENCE.
C
      DO 6 I=1,N
      T00=(G1(I)+G5(I))*DZOT
      T01(I)=(T00+DZ*G3(I))*.5
      T10(I)=(4.*T01(I)-T00)/3.
      IF (I.GT.3) GO TO 6
      TR=REAL(T01(I))
C      TI=AIMAG(T01(I))
      TI=IMAG(T01(I))
      TMAG1=TMAG1+TR*TR+TI*TI
      TR=REAL(T10(I))
C      TI=AIMAG(T10(I))
      TI=IMAG(T10(I))
      TMAG2=TMAG2+TR*TR+TI*TI
6     CONTINUE
      TMAG1=DSQRT(TMAG1)
      TMAG2=DSQRT(TMAG2)
      CALL TEST(TMAG1,TMAG2,TR,0.D0,0.D0,TI,DMIN)
      IF(TR.GT.RX)GO TO 8
      DO 7 I=1,N
7     SUM(I)=SUM(I)+T10(I)
      NT=NT+2
      GO TO 12
8     CALL SFLDS (Z+DZ*.25,G2)
      CALL SFLDS (Z+DZ*.75,G4)
      TMAG1=0.
      TMAG2=0.
C
C     EVALUATE 5 POINT ROMBERG RESULT AND TEST CONVERGENCE.
C
      DO 9 I=1,N
      T02=(T01(I)+DZOT*(G2(I)+G4(I)))*.5
      T11=(4.*T02-T01(I))/3.
      T20(I)=(16.*T11-T10(I))/15.
      IF (I.GT.3) GO TO 9
      TR=REAL(T11)
C      TI=AIMAG(T11)
      TI=IMAG(T11)
      TMAG1=TMAG1+TR*TR+TI*TI
      TR=REAL(T20(I))
C      TI=AIMAG(T20(I))
      TI=IMAG(T20(I))
      TMAG2=TMAG2+TR*TR+TI*TI
9     CONTINUE
      TMAG1=DSQRT(TMAG1)
      TMAG2=DSQRT(TMAG2)
      CALL TEST(TMAG1,TMAG2,TR,0.D0,0.D0,TI,DMIN)
      IF(TR.GT.RX)GO TO 14
10    DO 11 I=1,N
11    SUM(I)=SUM(I)+T20(I)
      NT=NT+1
12    Z=Z+DZ
      IF (Z.GT.ZEND) GO TO 17
      DO 13 I=1,N
13    G1(I)=G5(I)
      IF (NT.LT.NTS.OR.NS.LE.NX) GO TO 3
      NS=NS/2
      NT=1
      GO TO 3
14    NT=0
      IF (NS.LT.NM) GO TO 15
      WRITE(*,19)  Z
      GO TO 10
15    NS=NS*2
      DZ=S/NS
      DZOT=DZ*.5
      DO 16 I=1,N
      G5(I)=G3(I)
16    G3(I)=G2(I)
      GO TO 5
17    CONTINUE
C**
      RETURN
C
18    FORMAT (' ERROR - B LESS THAN A IN ROM2')
19    FORMAT (33H ROM2 -- STEP SIZE LIMITED AT Z =,1P,E12.5)
      END
C
C
C
      SUBROUTINE GX (ZZ,RH,XK,GZ,GZP)
C     SEGMENT END CONTRIBUTIONS FOR THIN WIRE APPROX.
      REAL*8 XK,R,R2,RK,CRK,SRK
      COMPLEX*16 GZ,GZP
C**
      R2=ZZ*ZZ+RH*RH
      R=DSQRT(R2)
      RK=XK*R
      CRK= DCOS(RK)
      SRK=-DSIN(RK)
      GZ=DCMPLX(CRK,SRK)/R
      GZP=-DCMPLX(1.,RK)*GZ/R2
      END
C
C
C
      SUBROUTINE GXX (ZZ,RH,A,A2,XK,IRA,G1,G1P,G2,G2P,G3,GZP)
C     SEGMENT END CONTRIBUTIONS FOR EXT. THIN WIRE APPROX.
      REAL*8 XK,R,R2,R4,RK,RK2,RH2
      COMPLEX*16 GZ,C1,C2,C3,G1,G1P,G2,G2P,G3,GZP
      R2=ZZ*ZZ+RH*RH
      R=DSQRT(R2)
      R4=R2*R2
      RK=XK*R
      RK2=RK*RK
      RH2=RH*RH
      T1=.25*A2*RH2/R4
      T2=.5*A2/R2
      C1=CMPLX(1.,RK)
      C2=3.*C1-RK2
      C3=CMPLX(6.,RK)*RK2-15.*C1
      GZ=CMPLX(DCOS(RK),-DSIN(RK))/R
      G2=GZ*(1.+T1*C2)
      G1=G2-T2*C1*GZ
      GZ=GZ/R2
      G2P=GZ*(T1*C3-C1)
      GZP=T2*C2*GZ
      G3=G2P+GZP
      G1P=G3*ZZ
      IF (IRA.EQ.1) GO TO 2
      G3=(G3+GZP)*RH
      GZP=-ZZ*C1*GZ
      IF (RH.GT.1.E-10) GO TO 1
      G2=0.
      G2P=0.
      RETURN
1     G2=G2/RH
      G2P=G2P*ZZ/RH
      RETURN
2     T2=.5*A
      G2=-T2*C1*GZ
      G2P=T2*GZ*C2/R2
      G3=RH2*G2P-A*GZ*C1
      G2P=G2P*ZZ
      GZP=-ZZ*C1*GZ
      RETURN
      END
C
C
C
      SUBROUTINE INTX (EL1,EL2,B,IJ,SGR,SGI)
C
C     INTX PERFORMS NUMERICAL INTEGRATION OF EXP(JKR)/R BY THE METHOD OF
C     VARIABLE INTERVAL WIDTH ROMBERG INTEGRATION.  THE INTEGRAND VALUE
C     IS SUPPLIED BY SUBROUTINE GF.
C
      INTEGER*4 NM,NS
      REAL*8 T01R,T10R,TE1R,T01I,T10I,TE1I,T11R,T20R,TE2R,TE2I,T11I,
     1 T20I,TE21,T00R,T00I,G1R,G1I,G2R,G2I,G3R,G3I,G4R,G4I,G5R,G5I,
     2 RKB2,B,S,EL1,EL2,Z,ZE,ZP,DZ,DZOT
      COMMON /TMI/ ZPK,RKB2,IJX
      DATA NX,NM,NTS,RX/1,65536,4,1.E-4/
C**
      Z=EL1
      ZE=EL2
      IF (IJ.EQ.0) ZE=0.
      S=ZE-Z
      FNM=NM
      EP=S/(10.*FNM)
      ZEND=ZE-EP
      SGR=0.
      SGI=0.
      NS=NX
      NT=0
C**
      CALL GF (Z,G1R,G1I)
C**
1     FNS=NS
      DZ=S/FNS
      ZP=Z+DZ
      IF (ZP-ZE) 3,3,2
2     DZ=ZE-Z
      IF (ABS(DZ)-EP) 17,17,3
3     DZOT=DZ*.5
      ZP=Z+DZOT
      CALL GF (ZP,G3R,G3I)
      ZP=Z+DZ
      CALL GF (ZP,G5R,G5I)
4     T00R=(G1R+G5R)*DZOT
      T00I=(G1I+G5I)*DZOT
      T01R=(T00R+DZ*G3R)*0.5
      T01I=(T00I+DZ*G3I)*0.5
      T10R=(4.0*T01R-T00R)/3.0
      T10I=(4.0*T01I-T00I)/3.0
C
C     TEST CONVERGENCE OF 3 POINT ROMBERG RESULT.
C
      CALL TEST (T01R,T10R,TE1R,T01I,T10I,TE1I,0.D0)
      IF (TE1I-RX) 5,5,6
5     IF (TE1R-RX) 8,8,6
6     ZP=Z+DZ*0.25
      CALL GF (ZP,G2R,G2I)
      ZP=Z+DZ*0.75
      CALL GF (ZP,G4R,G4I)
      T02R=(T01R+DZOT*(G2R+G4R))*0.5
      T02I=(T01I+DZOT*(G2I+G4I))*0.5
      T11R=(4.0*T02R-T01R)/3.0
      T11I=(4.0*T02I-T01I)/3.0
      T20R=(16.0*T11R-T10R)/15.0
      T20I=(16.0*T11I-T10I)/15.0
C
C     TEST CONVERGENCE OF 5 POINT ROMBERG RESULT.
C
      CALL TEST (T11R,T20R,TE2R,T11I,T20I,TE2I,0.D0)
      IF (TE2I-RX) 7,7,14
7     IF (TE2R-RX) 9,9,14
8     SGR=SGR+T10R
      SGI=SGI+T10I
      NT=NT+2
      GO TO 10
9     SGR=SGR+T20R
      SGI=SGI+T20I
      NT=NT+1
10    Z=Z+DZ
      IF (Z-ZEND) 11,17,17
11    G1R=G5R
      G1I=G5I
      IF (NT-NTS) 1,12,12
12    IF (NS-NX) 1,1,13
C
C     DOUBLE STEP SIZE
C
13    NS=NS/2
      NT=1
      GO TO 1
14    NT=0
      IF (NS-NM) 16,15,15
15    WRITE(*,20)  Z
      GO TO 9
C
C     HALVE STEP SIZE
C
16    NS=NS*2
      FNS=NS
      DZ=S/FNS
      DZOT=DZ*0.5
      G5R=G3R
      G5I=G3I
      G3R=G2R
      G3I=G2I
      GO TO 4
17    CONTINUE
      IF (IJ) 19,18,19
C
C     ADD CONTRIBUTION OF NEAR SINGULARITY FOR DIAGONAL TERM
C
18    SGR=2.*(SGR+DLOG((DSQRT(B*B+S*S)+S)/B))
      SGI=2.*SGI
19    CONTINUE
      RETURN
C
20    FORMAT (' STEP SIZE LIMITED AT Z=',F10.5)
      END
C
C
C
      SUBROUTINE GF (ZK,CO,SI)
C
C     GF COMPUTES THE INTEGRAND EXP(JKR)/(KR) FOR NUMERICAL INTEGRATION.
C
      REAL*8 ZK,CO,SI,RK,RKS,RKB2
      COMMON /TMI/ ZPK,RKB2,IJ
C**
      ZDK=ZK-ZPK
      RK=DSQRT(RKB2+ZDK*ZDK)
      SI=DSIN(RK)/RK
      IF (IJ) 1,2,1
1     CO=DCOS(RK)/RK
      RETURN
2     IF (RK.LT..2D0) GO TO 3
      CO=(DCOS(RK)-1.)/RK
      RETURN
3     RKS=RK*RK
      CO=((-1.38888889D-3*RKS+4.16666667D-2)*RKS-.5)*RK
      RETURN
      END
C
C
C
      SUBROUTINE TEST (F1R,F2R,TR,F1I,F2I,TI,DMIN)
C
C     TEST FOR CONVERGENCE IN NUMERICAL INTEGRATION
C
      REAL*8 F1R,F2R,TR,F1I,F2I,TI,DMIN,DEN
      DEN=ABS(F2R)
      TR=ABS(F2I)
      IF (DEN.LT.TR) DEN=TR
      IF (DEN.LT.DMIN) DEN=DMIN
      IF (DEN.LT.1.D-37) GO TO 1
      TR=ABS((F1R-F2R)/DEN)
      TI=ABS((F1I-F2I)/DEN)
      RETURN
1     TR=0.
      TI=0.
      RETURN
      END
C
C
C
      SUBROUTINE SFLDS (T,E)
C
C     SFLDX RETURNS THE FIELD DUE TO GROUND FOR A CURRENT ELEMENT ON
C     THE SOURCE SEGMENT AT T RELATIVE TO THE SEGMENT CENTER.
C
      REAL*8 T,PI,TP,THET,POT,RK,R1,R2,ZMH,ZPH,SFAC
      COMPLEX*16 E,U,U2,XX1,XX2,ERV,EZV,ERH,EZH,EPH
      COMPLEX*16 T1,EXK,EYK,EZK,EXS,EYS,EZS,EXC,EYC,EZC,
     1 ZRATI,ZRATI2,FRATI,ER,ET,HRV,HZV,HRH
      INTEGER*4 IND1,IND2
      COMMON/DATAJ/S,B,XJ,YJ,ZJ,CABJ,SABJ,SALPJ,EXK,EYK,EZK,EXS,EYS,
     1 EZS,EXC,EYC,EZC,RKH,IND1,IND2,IPGND,IEXK
      COMMON /INCOM/ XO,YO,ZO,SN,XSN,YSN,ISNOR
      COMMON /GWAV/ U,U2,XX1,XX2,R1,R2,ZMH,ZPH
      COMMON /GND/ZRATI,ZRATI2,FRATI,CL,CH,SCRWL,SCRWR,NRADL,KSYMP,
     1 IFAR,IPERF,T1,T2
      DIMENSION E(9)
      DATA PI/3.141592654D0/,TP/6.283185308D0/,POT/1.570796327D0/
C**
C     E      WRITE(*,*) '   SFLDS: START'
C**
      XT=XJ+T*CABJ
      YT=YJ+T*SABJ
      ZT=ZJ+T*SALPJ
      RHX=XO-XT
      RHY=YO-YT
      RHS=RHX*RHX+RHY*RHY
C      RHO=DSQRT(RHS)
      RHO=SQRT(RHS)
      IF (RHO.GT.0.) GO TO 1
      RHX=1.
      RHY=0.
      PHX=0.
      PHY=1.
      GO TO 2
1     RHX=RHX/RHO
      RHY=RHY/RHO
      PHX=-RHY
      PHY=RHX
2     CPH=RHX*XSN+RHY*YSN
      SPH=RHY*XSN-RHX*YSN
      IF (ABS(CPH).LT.1.E-10) CPH=0.
      IF (ABS(SPH).LT.1.E-10) SPH=0.
      ZPH=ZO+ZT
      ZPHS=ZPH*ZPH
      R2S=RHS+ZPHS
C      R2=DSQRT(R2S)
      R2=SQRT(R2S)
      RK=R2*TP
      XX2=CMPLX(DCOS(RK),-DSIN(RK))
      IF (ISNOR.EQ.1) GO TO 3
C
C     USE NORTON APPROXIMATION FOR FIELD DUE TO GROUND.  CURRENT IS
C     LUMPED AT SEGMENT CENTER WITH CURRENT MOMENT FOR CONSTANT, SINE,
C     OR COSINE DISTRIBUTION.
C
      ZMH=1.
      R1=1.
      XX1=0.
C**
      CALL GWAVE (ERV,EZV,ERH,EZH,EPH)
C**
      ET=-(0.,4.77134)*FRATI*XX2/(R2S*R2)
      ER=2.*ET*CMPLX(1.,RK)
      ET=ET*CMPLX(1.-RK*RK,RK)
      HRV=(ER+ET)*RHO*ZPH/R2S
      HZV=(ZPHS*ER-RHS*ET)/R2S
      HRH=(RHS*ER-ZPHS*ET)/R2S
      ERV=ERV-HRV
      EZV=EZV-HZV
      ERH=ERH+HRH
      EZH=EZH+HRV
      EPH=EPH+ET
      ERV=ERV*SALPJ
      EZV=EZV*SALPJ
      ERH=ERH*SN*CPH
      EZH=EZH*SN*CPH
      EPH=EPH*SN*SPH
      ERH=ERV+ERH
      E(1)=(ERH*RHX+EPH*PHX)*S
      E(2)=(ERH*RHY+EPH*PHY)*S
      E(3)=(EZV+EZH)*S
      E(4)=0.
      E(5)=0.
      E(6)=0.
      SFAC=PI*S
      SFAC=DSIN(SFAC)/SFAC
      E(7)=E(1)*SFAC
      E(8)=E(2)*SFAC
      E(9)=E(3)*SFAC
C**
C     E      WRITE(*,*) '   SFLDS: RETURN LINE 91'
C**
      RETURN
C
C     INTERPOLATE IN SOMMERFELD FIELD TABLES
C
3     IF (RHO.LT.1.E-12) GO TO 4
      THET=DATAN(ZPH/RHO)
      GO TO 5
4     THET=POT
5     CALL INTRP (R2,THET,ERV,EZV,ERH,EPH)
C     COMBINE VERTICAL AND HORIZONTAL COMPONENTS AND CONVERT TO X,Y,Z
C     COMPONENTS.  MULTIPLY BY EXP(-JKR)/R.
      XX2=XX2/R2
      SFAC=SN*CPH
      ERH=XX2*(SALPJ*ERV+SFAC*ERH)
      EZH=XX2*(SALPJ*EZV-SFAC*ERV)
      EPH=SN*SPH*XX2*EPH
C     X,Y,Z FIELDS FOR CONSTANT CURRENT
      E(1)=ERH*RHX+EPH*PHX
      E(2)=ERH*RHY+EPH*PHY
      E(3)=EZH
      RK=TP*T
C     X,Y,Z FIELDS FOR SINE CURRENT
      SFAC=DSIN(RK)
      E(4)=E(1)*SFAC
      E(5)=E(2)*SFAC
      E(6)=E(3)*SFAC
C     X,Y,Z FIELDS FOR COSINE CURRENT
      SFAC=DCOS(RK)
      E(7)=E(1)*SFAC
      E(8)=E(2)*SFAC
      E(9)=E(3)*SFAC
C**
C     E      WRITE(*,*) '   SFLDS: RETURN LINE 125'
C**
      RETURN
      END
